SoSe2022

Folienübersicht

Poweranalyse und Berechnung des Stichprobenumfangs

Übliche Experimentfragen

Nachdem das Design des Experiments und die statistische Analyse bestimmt wurden, ergeben sich oft folgende 3 Fragen:

  1. Welcher Stichprobenumfang ist nötig?
  2. Wie hoch ist die Teststärke oder Power des Experiments?
  3. Was ist der kleinste Effekt, der mit dem Experiment nachgewiesen werden kann?

Einflussgrößen

Für die Beantwortung dieser Fragen braucht es folgende 5 Informationen:

  1. Stichprobenumfang (N)
  2. Effektgröße (ES oder \(\delta\)) → die Mindestgröße, die nachweisbar ist oder die vorhergesagte Größe des Effekts.
    • Je nach Experiment handelt es sich um die Differenz zwischen 2 Mittelwerten (t-Test), Unterschiede in Proportionen, das Verhältnis zwischen ‘between-group’ und ‘within-group’ Variabilität (ANOVA) oder die Korrelation zwischen 2 Variablen.
  3. Variabilität des zu messenden Merkmals (\(\sigma^2\)) → bezieht sich beim t-Test auf die gepoolte Standardabweichung, bei der ANOVA auf die ‘within-group’ Variabilität.
  4. Teststärke (Power) → wenn diese fix sein soll, wird oft Wert von 0.8 oder 0.9 gewählt.
  5. Signifikanzniveau (\(\alpha\)) → gewöhnlich festgelegt auf 0.05.

Zusammenhang

Die universale Power-Gleichung

\[Power \propto \frac{ES~\alpha~\sqrt{n}}{\sigma^2}\]

Die Teststärke steigt proportional* mit

  • steigender Effektgröße
  • steigendem \(\alpha\) (wenn das Signifikanzniveau weniger strikt wird)
  • zunehmendem Stichprobenumfang
  • abnehmender Varianz

→ Durch Umstellung der Gleichung kann jede dieser 5 Größen berechnet werden, wenn die anderen 4 gegeben sind.

*Wie genau die Teststärke sich ändert hängt von der Analyse ab.

Zum Glück gibt es auch R Funktionen

  • Je nach statistischer Analyse gibt es eine spezifische Funktion mit der die Stichprobengröße (n), Teststärke (power) und testbare Effektgröße (delta) ermittelt werden kann, z.B:
    • t-Test: power.t.test()
    • Zweistichprobentest für Proportionen: power.prop.test()
    • 1-faktorielle ANOVA: power.anova.test()

Beispiel t-Test

power.t.test(n = NULL, delta = NULL, sd = 1, sig.level = 0.05,
             power = NULL,
             type = c("two.sample", "one.sample", "paired"),
             alternative = c("two.sided", "one.sided"),
             strict = FALSE, tol = .Machine$double.eps^0.25)

Berechnung des Stichprobenumfangs

t-Test

N beim unabhängigen t-Test

Um die Stichprobengröße für einen unabhängigen t-Test zu ermitteln, brauchen wir einfach nur die Formel der Teststatistik nach n umwandeln. Es gilt für gleiche Varianzen und Stichprobengrößen:

\[t=\frac{(\bar{X}_1-\bar{X}_2)}{\sqrt{2\frac{s_p^2}{n}}}=\frac{d}{\sqrt{2\frac{s_p^2}{n}}}~~\Rightarrow~~\sqrt{2\frac{s_p^2}{n}}=\frac{d}{t}~~\Rightarrow~~2\frac{s_p^2}{n}=\frac{d^2}{t^2}~~\Rightarrow~~n = \frac{2t^2s_p^2}{d^2}\]

  • \(s_p\) ist die gepoolte Standardabweichung bei gleicher Stichprobengröße: \(\sqrt{\frac{s_1^2+s_2^2}{2}}\)
  • \(t\) hängt von der Wahl der Power und des Signifikanzniveaus ab.
  • Werden diese auf 1-\(\beta\) = 0.8 und \(\alpha\) = 0.05 gesetzt, so beträgt t = 2.8 und \(t^2 = 7.84 \approx 8\) (die Quantilen einer Normalverteilung bei 0.05/2 und 0.8 werden addiert).
  • Somit gilt: \(n = \frac{2*8*s_p^2}{d^2}\)

Berechnung des Stichprobenumfangs

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Kenngröße Buchfink Mönchsgrasmücke
Mittelwert 1800km 3000km
Standardabweichung s ±900km ±1000km
Stichprobengröße n 20 30


Buchfink


Mönchsgrasmücke
d <- abs(1800-3000)
s_p <- sqrt((900^2 + 1000^2) / 2)
(n <- (16*s_p^2) / d^2)
[1] 10.05556
  • → Um einen Unterschied von 1200km statistisch nachweisen zu können, braucht es mindestens 10 Replikate pro Gruppe.

Berechnung des Stichprobenumfangs

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Kenngröße Buchfink Mönchsgrasmücke
Mittelwert 1800km 3000km
Standardabweichung s ±900km ±1000km
Stichprobengröße n 20 30


Buchfink


Mönchsgrasmücke
d <- abs(1800-3000)
s_p <- sqrt((900^2 + 1000^2) / 2)
(n <- (16*s_p^2) / d^2)
[1] 10.05556
  • → Um einen Unterschied von 1200km statistisch nachweisen zu können, braucht es mindestens 10 Replikate pro Gruppe.
power.t.test(power=0.8, sd=s_p, delta=d)
     Two-sample t test power calculation 

              n = 10.9148
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

Berechnung des Stichprobenumfangs

t-Test - Beispiel 2

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

Und wie groß müsste die Stichprobe sein, wenn wir bereits einen Unterschied von 600km bzw. 300km statistisch nachweisen wollen?

delta von 600km

power.t.test(power = 0.8, sd = s_p, delta = 600)
     Two-sample t test power calculation 

              n = 40.44558
          delta = 600
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

delta von 300km

power.t.test(power = 0.8, sd = s_p, delta = 300)
     Two-sample t test power calculation 

              n = 158.8158
          delta = 300
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

Shiny App zur Berechnung von N

Unabhängiger t-Test

Stichprobenkalkulation vs. Poweranalyse

Tuning der Stichprobengröße

  • Faktoren, die den Stichprobenumfang bestimmen:
    • Durchführbarkeit / Zeit / Konvention
  • Dadurch ist eigentlich oft schon im Vorfeld der Umfang festgelegt.
  • Aufgrund der oft überraschend hohen Stichprobengröße wird dann gerne so lange an den andere Stellschrauben gedreht, bis das erwünschte N raus kommt → das ist aber nicht Sinn und Zweck der Analyse.
  • Die Stichprobenanalyse wird daher am besten genutzt rauszufinden, ob das geplante Experiment Sinn macht oder nicht:
    • Wenn die Teststärke beim geplanten N nur noch bei 40% liegt, ist es unwahrscheinlich ein signifikantes Ergebnis zu erhalten - auch wenn es einen Effekt tatsächlich gibt.
  • Anstatt die Power festzulegen und N zu bestimmen, kann auch N auf einen maximal durchführbaren Wert festgelegt werden, um dann die Power zu berechnen.

Poweranalyse

t-Test - Beispiel 1

Beispiel Zugverhalten: Buchfinken vs. Mönchsgrasmücken

N = 20

power.t.test(n = 20, sd = s_p, delta = d)
     Two-sample t test power calculation 

              n = 20
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.9729651
    alternative = two.sided

NOTE: n is number in *each* group

Im vorliegenden Experiment betrug der Stichprobenumfang 20 und 30. Angenommen in beiden Fällen war N = 20, wie hoch ist dann die Teststärke gewesen?

N = 5

power.t.test(n = 5, sd = s_p, delta = d)
     Two-sample t test power calculation 

              n = 5
          delta = 1200
             sd = 951.3149
      sig.level = 0.05
          power = 0.4190796
    alternative = two.sided

NOTE: n is number in *each* group

Wie hoch wäre die Teststärke bei 5 Vogelmessungen gewesen, um einen Unterschied von 1200 km statistisch nachweisen zu können?

Shiny App zur Berechnung der Power

Unabhängiger t-Test

Berechnung des minimal nachweisbaren Effektes

t-Test - Beispiel

Angenommen es soll die Untersuchung des Zugverhaltens wiederholt werden, es können aber nur 8 Vögel jeweils untersucht werden. Welche Differenz könnte dann überhaupt statistisch nachgewiesen werden?

N = 5

power.t.test(power = 0.8, n = 8, sd = s_p)
     Two-sample t test power calculation 

              n = 8
          delta = 1433.285
             sd = 951.3149
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

→ Die Differenz müsste mindestens 1433.3 km betragen, um überhaupt einen signifikanten Unterschied zwischen den beiden Arten nachweisen zu können.

Sog. ‘Power Curves’

1

  • Anstatt immer nur einzelne Szenarios durchzuspielen ist es hilfreich, die Teststärke für mehrere Kombinationen an Parametern als Kurve zu visualisieren:
# Erstelle df mit versch. Kombinationen
df <- expand.grid(n = seq(5, 40, 2), sd = c(500,750,1000, 1250), 
  delta = c(500,1000), power = NA)
dim(df)
[1] 144   4
# Powerberechnung für jede der Kombinationen
for (i in 1:nrow(df)) {
  df$power[i] <- power.t.test(n = df$n[i], sd = df$sd[i], delta = df$delta[i])$power
}

Sog. ‘Power Curves’

2

df %>% ggplot(aes(n, power, group = sd, color = factor(sd))) +
  geom_line() + geom_hline(yintercept = 0.8) + facet_grid(~delta, labeller = label_both)

  • Bei hoher Effektgröße und geringer Variabilität wird die nicht-lineare Beziehung zwischen der Stichprobengröße und der Power sichtbar.
    • → Kleine Erhöhungen des Stichprobenumfangs führen hier bereits zu großen Erhöhungen der Teststärke.
  • Die Kurve macht schnell deutlich, dass ein N von 40 nicht ausreicht um eine Power von 0.8 zu erzielen, wenn die Standardabweichung >750 ist und die Differenz bei nur 500km liegt.

Simulations-basierte Poweranalyse

1

  • Die Effektgröße und die Variabilität sind 2 Parameter in der Poweranalyse, die meist recht ungewiss sind.
  • Daher sind die berechneten Stichprobengrößen und Powerschätzungen mit hohen Unsicherheiten belegt.
  • Anstelle zu fragen “Welchen Stichprobenumfang brauche ich um 80% Teststärke zu haben?” sollte besser gefragt werden “Welchen Stichprobenumfang brauche ich um 90% sicher zu sein, dass das Experiment eine Power von 80% hat?”

Simulation für das Zugverhalten-Beispiel

set.seed(123)
# Erzeugung von normalverteilten mittleren Differenzen und SDs
delta <- rnorm(5000, 1200, 100)
sigma <- rnorm(5000, 951, 100)
df <- data.frame(delta = delta, sd = sigma)

# Berechnung von N für jede Kombination an ES und SD
n <- apply(df, 1, function(x) power.t.test(delta = x[1], sd = x[2], power = 0.8)$n)

Simulations-basierte Poweranalyse

2

Wie hoch sollte N sein, um eine 90%ige Wahrscheinlichkeit zu haben, dass die Power 80% beträgt?

quantile(n, 0.9)
     90% 
14.91489 

→ Unsere ursprüngliche Poweranalyse empfahl eine N von 11.

p

Berechnung des Stichprobenumfangs

ANOVA

power.anova.test(groups = NULL, n = NULL,
                 between.var = NULL, within.var = NULL,
                 sig.level = 0.05, power = NULL)

Input-Paramater

  • Anzahl der Gruppen
  • Stichprobenumfang pro Gruppe (Annahme balanciertes Design)
  • Varianz zwischen den Gruppen
  • Varianz innerhalb der Gruppen
  • Signifikanzniveau
  • Teststärke

Berechnung des Stichprobenumfangs

ANOVA - Beispiel

Gewichtsunterschiede von Atlantischem Lachs, der in 4 verschiedenen Typen von Netzkäfigen gezüchtet wurde (N = 24).

df <- salmon %>%  group_by(cages) %>% 
  summarise(mean = mean(weight), 
    sd = sd(weight))
power.anova.test(groups = 4, 
  between.var = var(df$mean), 
  within.var = mean(df$sd)^2, 
  power = 0.8)
     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 2.345194
    between.var = 0.7663649
     within.var = 0.2335063
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group

Power-Kurve

ANOVA - Beispiel

n <- 2:10
p <- power.anova.test(groups = 4, n = n, 
  between.var = var(df$mean), 
  within.var = mean(df$sd)^2)
p
     Balanced one-way analysis of variance power calculation 

         groups = 4
              n = 2, 3, 4, 5, 6, 7, 8, 9, 10
    between.var = 0.7663649
     within.var = 0.2335063
      sig.level = 0.05
          power = 0.6193382, 0.9544039, 0.9968026, 0.9998332, 0.9999929, 0.9999997, 1.0000000, 1.0000000, 1.0000000

NOTE: n is number in each group
plot(n,p$power)

Berücksichtigungen

  • Die Poweranalyse und Stichprobenkalkulation basieren auf den gleichen Annahmen wie die zugrunde liegenden Tests:
    • Für eine zuverlässige Poweranalyse und Stichprobenkalkulation beim t-Test muss die Normalitätsannahme für jede Gruppe erfüllt sein.
    • Bei der Poweranalyse für eine ANOVA müssen alle Gruppen dazu die gleiche gemeinsame Varianz haben, damit dem Ergebnis vertraut werden kann.
  • Da die Effektgröße und Variabilität oft unsicher sind, nutze die Simulations-basierte Analyse.

Hilfreiches Paket: pwr

  • Funktionen für weitere Tests und Design gibt es in dem R Paket pwr von Stéphane Champely, welches sich an Cohen (1988) orientiert (siehe auch die Vignette):
Funktion Poweranalyse für
pwr.2p.test() 2 Proportionen (gleiches N)
pwr.2p2n.test() 2 Proportionen (ungleiches N)
pwr.anova.test() Balanced one-way ANOVA (gleiches N)
pwr.chisq.test() Chi-Quadrat-Test
pwr.f2.test() Generalisierte Lineare Modelle (GLM)
pwr.p.test() Proportion (1 Stichprobe)
pwr.r.test() Korrelation
pwr.t.test() t-Tests (1-2-Stichproben, unabhängig, gepaart)
pwr.t2n.test() t-Test (2-Stichproben, ungleiches N)

Weitere hilfreiche Pakete

  • powerAnalysis: Funktionen für einfache Tests (ähnlich wie pwr)
  • webpower: Sammlung von Funktionen für simple Tests bis komplexe lineare Modelle, SEMs, etc.
  • Superpower: Poweranalysen für komplexere ANOVA Designs
  • simglm: Spezialpaket für Poweranalysen mittels Simulationen
  • SIMR, SIMR, SIMR: Poweranalyse für Generalized Linear Mixed Effects Models

Übungsaufgabe

Übungen aus…

Kapitel 11 - Poweranalyse und Bestimmung des Stichprobenumfangs


  • R Notebook-Skripte
    • DS2_11_Übungen.Rmd
    • DS2_11_Übungen_Lösung.Rmd

Abschlussquiz

Fragen?